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The stability of dynamical states characterized by a uniform firing rate (splay states) is analyzed 
in a network of N globally pulse-coupled rotators (neurons) subject to a generic velocity field. 
In particular, we analyse short-wavelength modes that were known to be marginally stable in the 
infinite N limit and show that the corresponding Floquet exponent scale as 1/N 2 . Moreover, we find 
that the sign, and thereby the stability, of this spectral component is determined by the sign of the 
average derivative of the velocity field. For leaky-integrate-and-fire neurons, an analytic expression 
■ for the whole spectrum is obtained. In the intermediate case of continuous velocity fields, the Floquet 

exponents scale faster than 1/N 2 (namely, as 1/N 4 ) and we even find strictly neutral directions in 
a wider class than the sinusoidal velocity fields considered by Watanabe and Strogatz in Physica D 
(SJ ■ 74 (1994) 197-253. 
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I. INTRODUCTION 



Understanding the dynamical behaviour of highly interconnected systems is of primary importance for neural 
dynamics [![, metabolic cycles Q, cold atoms and synchronization in general oscillators [J. A wide variety of 
interesting phenomena has been discovered, but a detailed understanding is often lacking, to the extent that even the 
^ | stability properties of stationary states in globally coupled oscillators have not been fully clarified. 
j_J • In this paper we study an ensemble of iV identical rotators, i.e. dynamical systems characterized by a single 
dynamical variable, the "phase" x. This includes neural models of leaky-integrate-and-fire (LIF) type, since the 
variable x (the membrane potential) can be interpreted as a phase. This is done by identifying the maximum value of 
the potential (the spiking threshold, that, without loss of generality, we assume to be equal to 1) with the minimum 
(the resetting value assumed to be equal to - see the next section for further details), as if they corresponded to 
the angles 27r and 0. More precisely, we investigate the stability of splay states [H]. In a splay state all rotators 
[ follow the same periodic dynamics x(t) (x(t + T) = x(t)) but different time shifts that are evenly distributed (and 
1—1 1 take all multiples of T/N, modulus T). Splay states have been observed experimentally in multimode laser systems 
' 3| and electronic circuits Q. Numerical and theoretical analyses have been performed in Josephson junction arrays 
aj, globally coupled Ginzburg-Landau equations Q, globally coupled laser models Q, and pulse-coupled neuronal 
' networks [lOfl . In the context of neuronal networks, splay states have been also recently investigated in systems with 
S^Q , dynamic synapses [Tl| and in realistic neuronal models [f2j ■ 

The first detailed stability analysis of LIF neurons was performed by developing a mean-field approach that is based 
on the introduction of the probability distribution of the neuron phases ITol. Il3j | . The method is expected to work in 
£C) • infinite systems. More recently, another approach has been implemented [14( , which is based on the linearization of a 
suitable Poincare map and works for any number of oscillators. As a result, it has been discovered that the spectrum of 
0^ . Floquet exponents is composed of two components: (i) the growth rate of "long-wavelength" perturbations - perfectly 
identified also with the method described in Ref. [I(J] ; (ii) the growth rate of "short-wavelength" (SW) perturbations 
J> , that cannot be characterized with methods that involve a coarse-graining over small scales. As discussed in flij . 
the latter component plays a crucial role when the width of the transmitted pulses is comparable to or smaller than 
T/N, since it may give rise to instabilities of otherwise stable patterns. The same analysis has also revealed that for 
finite pulse-widths, SW are marginally stable in the infinite TV-limit. It is therefore important to investigate more 
thoroughly finite systems, because it is still unclear whether and when they are stable. 
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Here we address precisely this question, by first implementing a perturbative technique in the standard LIF model 
and by then numerically investigating the behaviour of a more general class of rotators, characterized by a nonlinear 
velocity field F(x) = x. All of our results indicate that the SW component scales as 1/N 2 if and only if F(l) F(Q). 
Moreover, we systematically find that the SW component is stable (resp. unstable) if F(l) < F(0) (resp. F(l) > F(0). 
Since AF = F(l) — F(0) is, by definition, the average derivative of F, the two classes of systems will be identified as 
decreasing and increasing fields, respectively. 

At the boundary between these two classes of fields, continuous velocity fields (F(l) = F(0)) turn out to exhibit a 
faster scaling to zero of the Floquet SW spectrum. In the case of analytic functions, many exponents appear even to 
be numerically indistinguishable from zero. This scenario is coherent with, and in some sense extends, the theorem 
proved in [l5l |. where it has been shown that in the presence of a sinusoidal field F — a(t) + sin(2irx + a), one should 
expect N — 3 zero exponents for any dependence of a(t). 

The paper is organized as follows. In section [TT] we introduce the model and the event-driven map that is used to 
carry out the stability analysis. In Sec. IIIII we derive analytical perturbative expressions for the Floquet spectrum in 
the case of LIF neurons. The results are compared with the numerical solution of the exact equation. In Sec. IIVI we 
numerically analyse several examples of velocity fields to test the validity of the above mentioned conjectures. Finally, 
in Sec. |Vl we summarize the main results and the open problems. 

II. THE MODEL 

We consider a network of N identical neurons (rotators) coupled via a mean-field term. The dynamics of the i-th 
neuron writes as 

±i = F{ Xi ) + gE{t) (1) 

where Xi represents the membrane potential, E{t) is the "mean" forcing field, and g is the coupling constant, the 
analysis will be limited to the excitatory case, i.e. g > 0. When the membrane potential reaches the threshold value 
Xi = 1, a spike is sent to all neurons (see below for the connection between single spikes and the global forcing field 
E) and it is reset to x% = 0. The resetting procedure is an approximate way to describe the discharge mechanism 
operating in real neurons. The function F{x) is assumed to be everywhere positive (thus ensuring that the neuron is 
repetitively firing, i.e. it is supra-threshold). For F(x) = a — x, the model reduces to the well known case of leaky 
integrate-and-fire (LIF) neurons. The field E is the linear superposition of the pulses emitted in the past when the 
membrane potential of each single neuron has reached the threshold value. By following Ref. [13], we assume that the 
shape of a pulse emitted at time t = is given by E s (t) = 2L±e~ at , where 1/a is the pulse-width. This is equivalent 
to saying that the total field evolves according to the equation 

2 

E(t) + 2aE(t) + a 2 E(t) = ^ £ 5{t - t n ) . (2) 

n\t n <t 

where the sum in the r.h.s. represents the source term due to the spikes emitted at times t n < t. 

A. Event-driven map 

As anticipated in the introduction, it is convenient to transform the differential equations into a discrete-time 
mapping. We do so by integrating Eq. ((21) from time t n to time t n +i (where t n is the time immediately after the n-th 
pulse has been emitted). The resulting map reads 

E(n + 1) = E(n)e- aT ^ + Q(n)r(n)e- aT ^ (3a) 

Q(n+l)=Q(n)e-^(n) + ?L , (3b) 

where r(n) = t n+ \ — t n is the interspike time interval and, for the sake of simplicity, we have introduced the new 
variable Q := aE + E . 

Moreover, the differential equation ((T|) can be formally integrated to obtain, 

x j (n+l)=F(x j (n),E{n),Q(n),T(n)) j = l,...,N —1 ; x m {n + 1) = (4) 



where m indicates the closest-to-threshold neuron at time n and the time interval r(n) is determined by imposing 
the condition that x m reaches the value 1 at time n + 1, immediately before being reset to zero. Altogether, we have 
therefore transformed the initial problem into a discrete time map for N + 1 variables: E, Q and N — 1 membrane 
potentials (where one degree of freedom is eliminated as a result of taking the Poincare-section) . A relevant property 
of identical mean-field coupled rotators is that the ordering of the local variables is preserved by the dynamical 
evolution: all neurons "rotate" around the circle [0,1] (1 being identified with 0) without passing each other. On 
the other hand, being the neurons identical, we can change their labels as they are indistinguishable. By following 
Ref. [3, [HI , it is convenient to start ordering the membrane potentials from the largest to the smallest one and then 
to introduce a comoving reference frame, i.e. to decrease by 1 the label of each neuron (plus 1 — > N) at each step of 
the iteration. In this frame, the label of the closest-to-threshold neuron is always equal to 1 and the splay state is just 
a fixed point of the transformation. Accordingly the linear stability analysis amounts to determining the eigenvalues 
of the corresponding linearized transformation. 

In order to carry out the stability analysis, it is necessary to derive an explicit expression for 
!F(xj(n), E(n), Q(n), t(ji)). This is not generally doable, but in the thermodynamic limit (N — ► oo) one can ex- 
ploit the smallness of t ~ 0(1/ N) and correspondingly set up a suitable perturbative expansion. We shall see that 
in order to correctly reproduce the stability of the splay state in typical cases, it is necessary to expand the map 
to fourth order. In a few peculiar models, a fully analytic calculation is possible. This is the case of LIF neurons, 
because of the linear structure of the velocity field: they will be analysed in the next section. 



III. LEAKY INTEGRATE-AND-FIRE MODEL 



Let us now consider the leaky integrate-and- fire case for the supra-threshold neuron, namely a > 1. In the comoving 
frame, Eq. (gj) writes as, [3, [lj] 

afj_i(n+l) =x j (n)e' T{n) + l-x 1 (n)e- T{n) j = 1, . . . , N - 1 , (5) 
with the boundary condition xn — 0, while the n-th intcrspikc interval is given by the self-consistent equation, 

a — xi(n) 



r(n) = In 
where, 



a + gH(n) - 1 



(6) 



H(n) = E(n) + 2±± - U Q( n ) . 7 

a— 1 \ a — 1 J (a — 1) 

In the absence of coupling (g — 0), e r is obviously equal to the ratio of the initial (a — x\) and final (a — 1) velocity 
of the first neuron, since the dynamics reduces to a pure relaxation. In that case there would be no way to determine 
r as the r.h.s. would be independent of it, but with interactions this is no longer true. As soon as the coupling is 
switched on, the velocity starts depending on the evolution of the field E in the way that it is summarized by the 
expression H(n). 

The set of Eqs. (|3I5I6I7P defines a discrete-time mapping that is perfectly equivalent to the original set of ordinary 
differential equations. It should be noticed that Eq. ([7]) is valid for any physically meaningful pulse-width value (i.e., 
a > 0) including a = 1, when there is no divergence or discontinuity. Moreover, in the parameter region considered 
in this paper (i.e. g > and a > 1) the logarithm in Eq. © is well defined, since one can show that H(n) is always 
positive. 

In this framework, the splay state reduces to a fixed point that satisfies the following conditions, 

r(n) = ^ , (8a) 

E(n) ee E , Q(n) = Q , (8b) 
Xj-! = x 3 e- T ' N + 1 - Xl e- T ' N , (8c) 

where T is the time elapsed between two consecutive spike emissions of the same neuron. A simple calculation yields, 

Q = ^ (l - e- aT ' N ) ~\e = tQ -I)'' . (9) 



The solution of Eq. (JSJi) involves a geometric series that, together with the boundary condition xn = , leads to a 
transcendental equation for the period T . This, in the large N limit and at the leading order, reduces to the following 
simple expression: 



T = In 



e- 1 - 
aT 



1 

- 9 



Xa-l)T + g\ 



(10a) 
(10b) 



The lack of any dependence of the period from the pulse-width is due to the fact that in the N — > oo limit the forcing 
field reduces to E = 1 /T [Kj, LLJ] . If we assume that a > 1 (which corresponds to assuming that the single neuron is 
supra-threshold) , we see that in the excitatory case (g > 0) the period T is well defined only for g < 1 (T — > 0, when 
g approaches 1), while in the inhibitory case (g < 0), a meaningful solution exists for any coupling strength (T — > oo 
for g — > — oo). 



A. Linear stability 

By linearizing Eqs. around the fixed point {5]), we obtain 

6E(n + 1) = e~ ar «;(7i) + Te~ aT 8Q(n) - (aE - Qe" Qr ) £r(n) , (11a) 
£Q(n + 1) = e- aT \SQ{n) - aQ<5r(n)j , (lib) 
5xj-i(n + 1) = e~ T [Sxj (n) — Sxi(n)] + e~ T (x\ — Xj)5r(n) , ( 11C ) 
and the expression for <5r(n) can be derived by linearizing Eqs. (|6|7p 

<5r(n) = r x fei (n) + T E SE{n) + T Q SQ(n) , (12) 

where t x :— dr / dx\ and analogous definitions are adopted for te and tq. 

In the comoving frame, the boundary condition xn = implies 5xn — 0. In practice, the stability problem is 
solved by computing the Floquet spectrum of multipliers {fJ-k} > k = 1, . . . , N + 1 corresponding to the linear evolution 
(fTTj) . It should be stressed that in general, the solution can be determined only numerically. 

However, it is convenient to rewrite the Floquet multipliers as 

( x k = e i ^e T ^ +iu ^ N , (13) 

where ipk — ^ff- , k = 1, . . . ,N — 1 and (fiN — ^Pn+i = 0, while Xk and uik are the real and imaginary parts of the 
Floquet exponents. The variable ipk plays the role of the wave- number in the linear stability analysis of spatially 
extended systems and one can say that Afc characterizes the stability of the fc-th mode. Previous studies [lj] have 
shown that the spectrum can be decomposed into two components depending on the index k: (i) k ~ 0(1); (ii) 
k/N ~ C(l). The first component corresponds to long- wavelength perturbations that can be formally analysed by 
taking the continuum limit (this was implicitely done in Ref. [10(] ) ; the second component corresponds to "high" 
frequency oscillations that require taking in full account the discreteness of the "spatial" index j. This is clearly 
illustrated in Fig. [1] where we have plotted the spatial component Sxj of (the real part of) three eigenvectors. The 
vector plotted in panel a) corresponds to (fk = 0.067T and is indeed both rather smooth and close to a sinusoidal 
function. In the other two panels, we can see that upon increasing the wave-number ipk, the discontinuous structure 
of the eigenvectors becomes increasingly evident. 

While the eigenvalues of the first component are of order 1 , the analysis carried out in Ref. [13] has revealed that 
the second component vanishes in the N — > oo limit. Therefore it is necessary to go beyond the zeroth order result 
to determine the stability of a splay state in large but finite system. As an example, in Fig. [5] we show the spectrum 
of the Floquet multipliers of the splay state for excitatory coupling (g > 0) and finite values of N . 



B. Analytical Results 

In the LIF model, many steps towards the determination of the Floquet exponents can be performed exactly. We 
start by deriving expressions that are valid for any number N of neurons and eventually introduce a perturbative 
approach to obtain an explicit expression in the large N limit. 
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FIG. 1: Three instances of the real part of eigenvectors for LIF neurons for a = 3, g = 0.4, a = 3, and N = 200. From top to 
bottom, panels (a-c) correspond to ip = 0.067T, 0.34-7T, and 0.787T. 
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FIG. 2: (Color online) Floquet exponent spectra for the LIF neurons: exact expression (|22|l (filled red squares), perturbative 
expression (|29p (blue line) and event driven map correct up to the fourth order in r (empty black circles). The parameters are 
a = 3.0, g = 0.4, and a = 30.0, N = 200. 



1. Exact Expressions 



We start by introducing the standard Ansatz 

SE(n + 1) = ii k 5E(n) 

From Eq. (fTTKb) 

aQ 



SQ 
6E 



V k e aT - 1 
arQ 



St 



{{i k e^ - l) 2 

By combining the above equations with Eq. (|12|). we find 



x T = £\ — a + 



gQe 



-(a-l)? 



5Q(n + 1) = fi k 5Q(n) 



aE — Q + aQr 
H k e aT - 1 



St 



a-l {a - l)(Ai fe e aT - 1) 



e, , g , ^ (1 - Mfee T ) 

a£; + ^T + ar0 (^^T) 



(14) 

(15) 
(16) 



(17) 



where x T denotes the derivative of x\ with respect to t. From the evolution equation for 8xj, Eq. (jllb). and by 
assuming that Sxj(n + 1) = /j,kSxj(n), we obtain 

fikfixj-i — e~ T (6xj — Sx\) + (x\ — Xj)e~ T 5t . (18) 

By using St — bx\jx T and introducing the expression (|5fc) for ij, we find 

e O~i)T § Xl / i \ 
^-1=6-^ + ^-^— -^- w - Y)+e -^6x 1 . (19) 

The solution of this recursive equation reads 



5 Xj = - ( I + e-A 5X1 + u + K & T • (20) 



We can determine the constant K by imposing that the above equation is an identity for j = 1. As a result, 

ft - (-1 + .-r) 1 - , f , ■ - ' + <-e"-'" (2D 

fe x \x T (e T -l) J Mfe-e" r a; r (e T -l) /i fe - 1 ^ k 

The equation for the determinant is finally obtained by imposing Sxp? = 0, 



r-T _ N-l I _ W-l 

* r (e T - l)^" 1 = - (z T (e T - 1) + e') ^- ^- + (22) 



Eq. (|22[l is an exact but implicit expression for all Floquet multipliers that applies for a generic number N of neurons. 
A numerical solution of Eq. ([22|) reveals that, for finite N, excitatory coupling, and a smaller than the critical value 
«c = dc{g,N) (see alsofijl), the splay state is strictly stable, although the maximum Floquet exponent approaches 
zero for increasing N [14j. In fact, as shown in Ref. pTj, in the limit N — ► oo, SW modes are marginally stable, i.e. 
X k = u>k = 0. 



2. Perturbative Expansion 

Since the Floquet exponents of the SW vectors are exactly equal to zero in the infinite N- limit, it is natural to 
investigate the stability of finite systems in a perturbative way. In particular, we find it convenient to introduce the 
smallness parameter r ~ 1/N. A posteriori, and in agreement with the numerical observations in [14(, it turns out 
that an expansion up to second order in r is necessary and sufficient to correctly determine the leading contribution 
of the Floquet spectrum. 

Let us start by expanding the stationary solution for Q and E ^ 



Q 



a 
T 



1 



E 



1 

T 



a 



a 
12 



12 



(23a) 
(23b) 



Next, we express the Floquet multipliers as 



(24) 



which amounts to assume that the Floquet exponent is proportional to r 2 . 
By expanding x T and with the help of Eq. (|10b ). we obtain 



1 + t - Bt 1 



(25) 



where 



1 qo? 
B = — + - — 

2 T 



1 

12 



1) 



(g«Vfc — I) 2 



(26) 



After inserting the above expansions into Eq. (|22p , we obtain 



(-l-T + Br*)M-*-tf) = -T*[B 



= ) + e-(l-/xf- 1 ) 



Now, with the help of Eq. ([24)) and completing the r expansion, 



1 



r fc T = ( B + - ) (1 - e 



(27) 



(28) 



By finally replacing the definition 
spectrum, 



of B into the above equation, we obtain an explicit expression of the Floquet 



Afe 

t2 



12T 2 



6 



(cOS^fc - 1) 



(29) 



In Fig. [21 Eq. ([290 is compared with the numerical but exact solution of Eq. ([2"2"| for N = 200, revealing a very good 
agreement. The divergence of this expression for a — > oo indicates that the SW modes are characterized by a different 
scaling behaviour for J- like pulses. In fact, as shown in Ref. [HI], the corresponding exponents do not scale to zero 
for N — > oo. Moreover, it is worth recalling that the stability of networks with strictly <5-pulses cannot be inferred by 
taking the limit a — > oo of the exact (non-perturbative) expressions [lij . 



IV. GENERAL CASE 



In the previous section we have seen that in the LIF model the SW component of the spectrum scales as 1/iV 2 and 
obtained an analytic expression for the leading term. It is natural to ask whether the observed scaling behaviour is 
peculiar to this system or it is a g eneral characteristics of pulse coupled oscillators. For F(x) = a + sin(27ra; + a), 
the general theorem proved in [lj| tells us that the dynamics of the oscillators is characterized precisely by N — 3 
zero exponents independently of the behaviour of the forcing field E. Therefore, it is an example of perfect neutral 
stability for any N. For generic velocity fields, it is not possible to obtain analytic expression, so that one has to rely 
on approximate expressions. Since we expect Afe — > for increasing N, it is natural to follow a perturbative approach 
from the very beginning, i.e. from the definition of the event-driven map. In [l4j], it has been shown that a second 
order expansion fails to reproduce the stability properties of the LIF model even on a qualitative level. In fact, the 
Floquet exponent of the single-step map is of order t 3 . This would naively suggest that a third order is sufficient; 
however, the eigenvalue equation involves an "integration" over N steps. Therefore, it is necessary to control the 
accuracy of the single iterate of the map up to order 1/-/V 4 , what is incidentally guaranteed by standard integration 
algorithms like fourth-order Runge-Kutta. Here he have preferred to determine an explicit expression for the map to 
be thereby linearized and used to determine the entire Floquet spectrum. The results plotted in Fig. [2] confirm that 
upon including terms up to 0(l/iV 4 ), we are able to reproduce the expected results. 

With the goal of identifying the typical scaling behaviour of the Floquet spectrum, in the following we investigate 
various types of functions F(x), starting from smooth periodic functions. As mentioned above, if F{x) contains just 



i.e. 

a 



if F(x) = a — sin(27ra;) , we expect N — 3 zero exponents 

30 and different numbers of rotators. The data indeed show that 4 exponents remain finite 



15J . In Fig. [3k., we plot the spectrum 



the main harmonic, 
for a = 3 , g = 0.4, 

for N — > oo, in agreement with the the theoretical results [Hj], since in our system there are N+l degrees of freedom. 
Moreover, we can see that the vast majority of the exponents are equal to zero within numerical accuracy. This is 
even beyond our expectations, because of the finite (fourth order) accuracy of the numerical computations. 

In the presence of more harmonics, there are no theoretical predictions which can guide us. In Fig. [3Jd, we present 
the results for F(x) = 3 — sin(27rx)/2 — 0.1sin(47rx) — 0.01 sin(67rx). There we see that there are no substantial 
differences from the previous case, the main novelty being that the number of (negative) exponents which remain 
finite for N — > oo is definitely larger than 4, probably 24 or 26 with our numerical accuracy. The presence of many 
zero exponents is confirmed by simulations performed for different parameter values. Whether the "numerical zeros" 
correspond to exact zeros and thereby to some conservation laws is however an open question. 

The choice of a periodic function F(x) such as a — sin(27ra;) is natural in the context of coupled rotators, where x 
is a true phase and the 0, 1 values can be identified with one another as they correspond to angles differing by 27r. 
In the LIF model, and I correspond to two different membrane potentials (actually, the minimum and maximum 
accessible values) and there is no reason a priori to expect F(l) = F(0), namely AF = F(l) — F(0) = —I. It is 
therefore important to understand whether the different scaling behaviour is to be attributed to the presence of a 
"discontinuity" in the velocity field and if the sign of the difference F(l) — F(0) matters or not. In order to clarify this 
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FIG. 3: (Color online) Floquet spectra: (a) single harmonic; (b) three harmonics. The parameters are a = 3.0, g = 0.4, and 
a = 30.0. 




FIG. 4: (Color online) (a) The two considered velocity fields: upper and lower curves correspond to: F(x) = 1.3 + 0.7x — x 2 
(solid line) and F(x) — 1.3 — 1.3a; + x 2 (dashed line); (b) the lower (resp. upper) Floquet spectra are associated to the lower 
(resp. upper) velocity field in Fig. |3Ji, please notice that in the inset due to the negative values of the spectra the correspondence 
is reversed. The values of the other parameters are g = 0.4, and a = 6.0. 



point we have investigated two parabolic fields with opposite concavities, but identical (and negative) nonzero value 
of the velocity difference at the extrema of the definition interval, i.e. AF = —0.3 (see Fig. H^, for their graphical 
representation). The results plotted in Fig. [JJd indicate that the presence of nonlinearities in the velocity field do 
neither affect the scaling of the spectrum, that is still proportional to r 2 , nor the overall stability: the whole branch 
is strictly negative. 

As a next step, we have investigated two increasing parabolic velocity fields with opposite concavities, but identical 
and positive difference at the extrema AF = 0.3 (see Fig. [SJi for a graphical representation). The results plotted in 
Fig. ]Ejp confirm once again the t 2 scaling. However, the stability has changed: now the SW component is positive, 
indicating that the splay state is weakly unstable. Altogether, we can summarize the results under the conjecture that 
all discontinuous velocity fields (i.e., where F(l) ^ F(0)) are characterized by exponents that scale as t 2 . Moreover, 
the stability depends on whether on the average the field increases or decreases. The analysis of several other velocity 
fields has confirmed this conjecture. 

In between the two classes of increasing and decreasing fields, there are continuous fields. The neutral stability 
of the sinusoidal fields is logically consistent with the observation that the stability depends on the sign of AF. 
In order to further explore the generality of the scenario, we have analysed two other cases: (i) a parabolic field 
F(x) = 1.3 - x(x - 1); (ii) a sinusoidal field F(x) = 1.3 - 0.25sin(7ra). Both fields are C (0) , but not C^, since 
the derivative in x = 1 and x = differ from one another. In Fig. \E\ we see that the spectra scale as 1/N 4 . This 




FIG. 5: (Color online) (a) The two considered velocity fields: upper and lower curves correspond to: F(x) = 1.3 + 1.3x — x 2 
and F(x) = 1.3 — 0.7a; + x 2 ; (b) the lower (resp. upper) Floquet spectra are associated to the lower (resp. upper) velocity field 
in Fig. [5K- The values of the other parameters are g = 0.4, and a = 6.0. 
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FIG. 6: (Color online) Floquet spectra: (a) continuous parabolic field F(x) — 1.3 — x(x — 1); (b) continuous sinusoidal field 
F(x) — 1.3 — 0.25 sin(7ra;) with period 2. The data refer to g — 0.4 and a = 6 and have been obtained by employing an 
event-driven map including terms up to order O(oo/Af v ). 



confirms that continuous functions exhibit an intermediate behaviour between positive and negative discontinuities. 
The scaling behaviour, as r 4 , has been verified by employing an event-driven map correct to order r 5 . Moreover, it is 
also interesting to notice the difference with respect to the analytic sinusoidal functions. In fact, it seems that exactly 
zero exponents are detected only for analytic velocity fields. 

To further verify the scenario, we have introduced a velocity field 

F(w) =a + &w(w — l) 

w = (af +x)/2 (30) 

parametrized by the exponent 7. The function is periodic but, upon increasing 7, it becomes increasingly steeper in the 
vicinity of x — 1 as shown in Fig. [7]^ a). For 7 < 10 we observe the same behaviour found for other periodic functions, 
i.e. the spectrum scales rapidly to zero (like 1/iV 4 ). To exemplify this case let us consider 7 = 2, as shown in Fig. 03 
the eigenvalue Vm/2 decreases as ~ I/TV 2 for sufficiently large N . Please notice that the eigenvalues are approaching 
zero from positive values in this case. For 7 ~ 10 2 the situation becomes more complicated: for sufficiently small 
TV- values T N / 2 is negative and almost constant (indicating an l/N 2 scaling of the Floquet eigenvalues), while by 
increasing N it becomes positive and Tjv/2 - ¥ only for very large N, — > (see Fig. [Th). By increasing 7, 
the 1/N 2 scaling region widens, but the overall behaviour is maintained, as shown in Fig. [7ji for 7 = 10,000. This 
suggests that for not too large values of N, the system perceives the field as if it was discontinuous, while at large 
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FIG. 7: (Color online) (a) The velocity fields expressed by Eqs. (|30|l for three different 7 values are reported in proximity of 
x = 1. The stability exponent T N / 2 corresponding to the highest frequency (ifiN/2 = is plotted in panels (b), (c), and (d) 
for the three different 7 values as a function of N. In (b) , the (black) filled circles refer to 7 = 2 and the dashed (red) line 
represents a power-law decay N -13 with (3 = 2. In (c), [rjv/2] 1 ^ 3 is displayed for 7 = 100; the dashed (blue) line indicates the 
zero axis. Finally, for a better visualization of the data for 7 = 10, 000, the vertical scale in (d) has been obtained by first 
shifting Tjv/2 by 20 units and thereby using exponentially separated units. All eigenvalues refer to a — 1.3, g = 0.4 and a = 6. 



N it crosses to continuous fields. This crossing is joined to a change of sign from negative values (as expected for 
discontinous fields with AF < 0) to positive ones. 

Finally, for the sake of mathematical generality, we have investigated the role of an additional, intermediate, 
discontinuity in the velocity field. More precisely, we have examined the piece-wise linear model, 

f a - 6/2 - for x < 0.5 . 
[X) \ a + 6/2 - x, for x > 0.5 > 

with a discontinuity of size 6 at x = 0.5. As shown in Fig. [8^, as a result of the discontinuity, the Floquet spectrum 
widens to cover a thick band, that is filled in an increasingly uniform way, upon increasing N . Depending whether the 
discontinuity is positive or negative, the band develops towards higher or smaller values, respectively (see Fig. [SJd), 
while the standard LIF spectrum (solid line) represents the locus of minima or maxima, respectively of such bands. 
Due to the finite thickness of the band itself, unstable SW modes can appear even for AF < 0. Finally, for 6=1, 
when AF = F(l) — F(0) = 0, the SW modes scale faster than 1/N 2 as for standard continuous functions. Altogether, 
the presence of an additional discontinuity does not modify the scaling behaviour, that is still controlled by the sign 
of AF. 



V. CONCLUSIONS AND OPEN PROBLEMS 



In this paper we have shown that the stability of splay states in pulse-coupled oscillators with generic velocity fields 
F(x) can be determined by rewriting the dynamics as event-driven maps. In particular we focused our attention 
on the SW modes. For discontinuous velocity fields, like that associated to leaky integrate-and-fire neurons, we find 




FIG. 8: (Color online) Floquet spectra I\ associated to the velocity field (|3ip . (a) Spectra for various N with b = +0.2: 
N = 800 (red) circles and N = 1,600 (blue) crosses, (b) Spectra for N = 1,600 at different fr-values: namely, (red) squares 
refer to b = —0.2, (green) diamonds to b = +0.2 and (blue) stars to b = +0.6. For comparison also the Floquet spectrum 
corresponding to a LIF neuron with the same parameters is reported (black solid line). All the data refer to a = 1.7, g = 0.4 
and a = 6. 

that all SW modes are stable, when the field on the average decreases (i.e. AF < 0), and unstable in the opposite 
case. Notice that this weak instability cannot be captured by the mean-field approach introduced in [lolh as the 
coarse-graining washes out SW modes. It is instructive to compare the role of AF, with the results of Refs. [13, (3- 
While studying an excitatory network of globally coupled rotators (in the limit of 5-like pulses) , Mirollo and Strogatz 
proved that the synchronous state is fully stable when x(<j>) is concave-down [13], where the phase <j) is nothing but the 
time variable (apart from a scaling factor) and the evolution refers to the single neuron dynamics x = F{x). Recently, 
these results have been complemented by the observation that clustered states are stable (while the synchronous 
regime is unstable), when x((f>) is concave-up [l8L Il9j . It is easy to verify that in a linear LIF neuron, the change of 
concavity is one-to-one connected with a change of sign of AF. More in general, a concave-up (-down) x(t) implies 
that AF > (AF < 0), while the opposite implication does not hold. For instance, for F(x) = 1.3 — sin[7r(9x — l)/4] 
(x E [0.1]), x(t) exhibits even two changes of concavity and yet we find that the SW modes are stable and scale as 
1/N 2 , as expected, since AF = — \/2/2. Altogether, the condition arising from the sign of F is more general than 
that based on the sign of the concavity of x(t), but it refers to SW modes only. 

Naively, one might think that that our results follow from the fact that AF is the average derivative of the velocity 
field in the interval [0, 1]. If, on the average, x = (AF)x, it is natural to expect exponential instability when AF > 0, 
stability in the opposite case, and marginal stability for AF = 0. However, to leading order, all SW modes are 
"marginally" stable. Presumably, there is some truth in the argument, but some refinements are required to put 
it on a firm basis. Furthermore, the 1/N 2 scaling of the Floquet spectrum has been so far rigorously proved for 
LIF neurons and has been confirmed by the numerical analysis of several nonlinear fields. This includes the well 
known exponential integrate-and-fire neurons (EIF) [2l| which, having a discontinuous velocity field with AF < 0, 
are consistently characterized by stable SW modes. 

The intermediate case of continuous velocity fields presents even more subtleties. The are generically characterized 
by a 1/N 4 scaling to zero, but we are unable to conclude whether the scaling law is entirely determined by the 
analyticity properties of the velocity field. Anyway, for fields composed of a few harmonics, it appears that all 
Floquet exponents (with the exception of a finite number of them) are equal to zero. These results suggest that 
perfectly marginal modes exist in a wide range of cases than that those proved in [151 ] and shown in [221 ] . The question 
is not of purely academic interest, since analytic velocity fields are generically encountered when dealing with coupled 
rotators, where x is a true phase. Furthermore, the so-called quadratic integrate-and-fire neurons (QIF) [23| belong 
to this class, as the velocity field is F(x) = dx{l — x) + e and it is therefore continuous. We have verified that this 
model is not an exception, as its short-wavelength spectrum scales faster than 1/N 2 . 
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